% syms r q d q n k
% for i=1:2
%     mu(i,1)=sym(strcat('mu',int2str(i)),'real');
%     n(i,1)=sym(strcat('n',int2str(i)),'real');
%     k(i,1)=sym(strcat('k',int2str(i)),'real');
% end
% mu=k'*n
% r=(1+q)/mu
% Dr_n = jacobian(r, n);    % -1*k'*(q+1)/(k'*n)^2
% Dr_q = jacobian(r, q);           % 1/(k'*n)
% % Dq_n = -K
% % Dq_q =  1

%h=double(h);
% pretty(B)

% ------------------------------------
classdef  test_linearize
    properties
    end
    methods
         function obj=test_linearize()
         end
        
         function   [Dr_n, Dr_q, Dq_n, Dq_q,A,B] =  coefs(obj,k,q,n)
                    Dr_n = -1*k'*(q+1)/(k'*n)^2; 
                    Dr_q =  1/(k'*n); 
                    Dq_n = -1.*k'; 
                    Dq_q =  1; 
                    A= [Dq_q  0
                          Dr_q    0];
                     B=[Dq_n
                          Dr_n];  
         end
         
        function d_yp=model_out(obj,k,u,n,q, demand)  
                      k_ = k *  prices_amz().get_compute_units/demand;
                      [Dr_n, Dr_q, Dq_n, Dq_q,A,B]  = obj.coefs(k_,q,n);    
                  
                      mo = ss(Dq_q,Dq_n,Dr_q,Dr_n,1);
                     d_yp = lsim(mo,u);  
        end

        function [k_hat, q,n,m_rt,x0] = test(obj, m_inst_count_vec)
            
            obj_k = test_mpc2_kalman_integral();
            load('data/experiment-60min.mat'); 
            % set the instances so around the desired point, allocatet your
            % self , 2 sec desire
            
            users = 800; 
            [m_rt , m_thr] = obj_k.get_steady_state_output_for( m_inst_count_vec , users); 
            [  d_rt_vec]= obj_k.linearize_around(rt_vec , m_rt);
            [ d_inst_count_vec] = obj_k.linearize_around(inst_count_vec,...
                m_inst_count_vec);
            n= m_inst_count_vec';
            u =  d_inst_count_vec;
            y = d_rt_vec;
            q= m_rt * m_thr;
            demand = 0.09; 
            k0 = 1;
            % prices_amz().get_compute_units/demand; 
           % k0 = ones(11,1); 
            
            mdl = @(k,u)model_out(obj,k,u,m_inst_count_vec',q,demand)    
            [k_hat,r,J,cov,mse] = nlinfit(u,y,mdl,k0);
            
            
            %--------------------------
            [Dr_n, Dr_q, Dq_n, Dq_q,A,B] = obj.coefs(k_hat,q,m_inst_count_vec'); 
            C = [0 1];  D=zeros(1,11);
            data=iddata(y, u , 1);           
            m=idss(A,B,C,D,0,'Ts', 1,'DisturbanceModel','Estimate','InitialState','Estimate');  % None   
            x0 = findstates(m,data); 
        end

        %----------------------------------------
             function d_yp=model_out2(obj,k,u)                        
                    k_ = k(1) *  prices_amz().get_compute_units/demand;
                    lambda = k(2); 
                    x=[0;0]
                    for i=1:size(u,2)
                        x(i+1,1)=x(i,1)+lambda-k_*u(i);      %q
                        x(i+1,2)=x(i,1)/(k_)
                    end
                      
                  
                      mo = ss(Dq_q,Dq_n,Dr_q,Dr_n,1);
                     d_yp = lsim(mo,u);  
             end
        
             function y=mdl(obj,a,x)
                a
                loga=a;
             %   y = exp(loga(1)) + exp(loga(2))*exp(-exp(loga(3))*x);  
            %    y = a(1) + a(2)*exp(-a(3)*x);  
                y=1./(a(1)*x-a(2)); 
             end
             
           function [k_hat, q,n,m_rt,x0] = test2(obj)
            format('shortG');
            warning('OFF');
            obj_k = test_mpc2_kalman_integral();
            load('data/experiment-60min.mat'); 
            % set the instances so around the desired point, allocatet your
            % self , 2 sec desire
            
        %    users = 800; 
            rt=[]; thr=[];  q=[];  users=[];
            inst_count_vec=[]; 
%              inst_count_vec=[...
%                  1 0 0 0 0 0 0 0 0 0 0 
%                  2 0 0 0 0 0 0 0 0 0 0 
%                  2 0 0 0 0 0 0 0 0 0 0 
%                  3 0 0 0 0 0 0 0 0 0 0 ];
%             
%              
%              users(:,1)=[...
%                  800
%                  850
%                  910
%                  780];
            cmp=prices_amz().get_compute_units(); 
            k=cumsum(randn(200,1) ); k= k-min(k); k=k./max(k); k=k*3+1; k=round(k);
            cu = round(sin((1:50)/5)*20+21); 
            for i=1:10  
                users(i,1)=800+round(randn*50);
               index=mod(round(rand*4+0.5),4)+1; 
             %  index=k(i);
                inst_count_vec(i,:) = zeros(1,11); 
             %   inst_count_vec(i,index) = ceil(cu(i)./cmp(index)); 
                inst_count_vec(i,index) =round(rand*20+0.5);                
                [rt(i,1) , thr(i,1)] = obj_k.get_steady_state_output_for(inst_count_vec(i,:) , users(i,1));     % q= m_rt * m_thr;  %demand = 0.09; 
                 q(i,1)=rt(i,1) * thr(i,1); 
                 disp(sprintf('i=%d',i));
            end  
 
%            (z+r)* n*mu/(n+k-1) = n
%             (z+r)*mu=n+k-1
            
pu=inst_count_vec(1:end,1:4)*cmp(1:4);
subplot(2,1,1); scatter(pu, rt(1:end));   
% mdl = @(a,x)(a(1) + a(2)*exp(-a(3)*x));
cvx_begin
   variables a b
   a*pu-b > 0;        
%   minimize( norm(rtt-b(:,3)) );       
   minimize( norm(a*pu-b-(1./rt) ));       
 cvx_end
xrange = min(pu):1:max(pu);
hold on
plot(xrange,1./(a*xrange-b))
%  a =     0.016734
%  b =     0.025036

% a0 = [2;2;2];
% [ahat,r,J,cov,mse] = nlinfit(pu,rt,@obj.mdl,a0)
% xrange = min(pu):1:max(pu);
% hold on
% plot(xrange,obj.mdl(ahat,xrange),'r')
% hold off


subplot(2,1,2); scatter(pu(2:end)-pu(1:end-1), rt(2:end)- rt(1:end-1));
 %-------------------------------------------
 % inst_count_vec(2:end,1:4)
 A=[ inst_count_vec(2:end,1:4)*cmp(1:4) thr(1:end-1) rt(1:end-1)] %x1:cu(t) x2:rt(t)
 b=[thr(2:end) rt(2:end)]
 x1 = inv(A'*A)*A'*b
bv=A*x1;
 disp(sprintf('achieved norm: %d',    norm(b(:,1)-bv(:,1))));                       %mean(rt)
 [inst_count_vec(2:end,1:4) inst_count_vec(2:end,1:4)*cmp(1:4) thr(2:end) rt(2:end)]
 [ inst_count_vec(2:end,1:4) b(:,1) bv(:,1)]
  inst_count_vec(2:end,1:4)*cmp(1:4)
 
 x1=x1'; 
 A=x1(:,2:3);  B=x1(:,1);
 x0=[ 92.669    ;   9.0185]   %thr rt 
  mo = ss(A, B*cmp(1:4)', [0 1],0,1);
  d_yp=zeros(10,4);
  for i=1:4
      disp(sprintf('i=%d',i*4));
      u=ones(10,1)*[i*4 0 0 0];
      d_yp(:,i) = lsim(mo,u,(1:10)',x0);
  end
  d_yp
 
  %---------------------------------
 %   A = ones(        B=zeros(1,1);
    C = [1 0];  D=zeros(1,1);
    
    u=A; y=b; 
    data=iddata(y, u , 1);           
    m=idss(A,B,C,D,0,'Ts', 1,'DisturbanceModel','Estimate','InitialState','Estimate');  % None   
    x0 = findstates(m,data); 
  
  %----------------------------------
  n=size(b,1); 
cvx_begin
   variables x(5,1) bv(n,1) rtt(n,1);
   bv == A*x;
   rtt== bv(:,1);      
   A*x > 0;        
%   minimize( norm(rtt-b(:,3)) );       
   minimize( norm(A*x-b) );       
 cvx_end
x
 [b(:,3) bv(:,3)]

            inv(A'*A)*A'*b
                prices_amz().get_compute_units
           kk = m_inst_count_vec .* prices_amz().get_compute_units'
           (kk/sum(kk))*104
        %   m_thr  m_inst_count_vec*k
           
 
           
          u =  inst_count_vec;
            y = rt_vec;
        
            
            k0 = 1;
            % prices_amz().get_compute_units/demand; 
           % k0 = ones(11,1); 
            
%             mdl = @(k,u)model_out(obj,k,u,m_inst_count_vec',q,demand)    
%             [k_hat,r,J,cov,mse] = nlinfit(u,y,mdl,k0);
%             
            
            %--------------------------
            [Dr_n, Dr_q, Dq_n, Dq_q,A,B] = obj.coefs(k_hat,q,m_inst_count_vec'); 
            C = [0 1];  D=zeros(1,11);
            data=iddata(y, u , 1);           
            m=idss(A,B,C,D,0,'Ts', 1,'DisturbanceModel','Estimate','InitialState','Estimate');  % None   
            x0 = findstates(m,data); 

        end

        % % X = [del_q ; del_r] 
        % T = size(u',2); 
        % n=m_inst_count_vec';
        % cvx_begin 
        %      variables k(11,1) X(2,T+1) A(2,2) B(2,11) Dq_q(1,1) Dr_q(1,1)  Dq_n(1,11) Dr_n(1,11) 
        %    %  Dr_n == -1*k'*(q+1)/(k'*n)^2; 
        %  %     Dr_q ==  1/(k'*n); 
        % %      Dq_n == -1.*k'; 
        % %      Dq_q ==  1; 
        % %      A== [Dq_q  0
        % %                Dr_q    0];
        % %      B==[Dq_n
        % %               Dr_n]; 
        %      X(2,1:T) == y';
        %     % minimize( norm( X(:,2:T+1) - [1  0; 1/(k'*n)   0] *X(:,1:T)+ [-1.*k'  ;  -1*k'*(q+1)/(k'*n)^2] * u'))
        %      minimize( norm( X(:,2:T+1) - [1  0; quad_over_lin(1,k'*n)   0] *X(:,1:T)+ [-1.*k'  ;  -1*k'*(q+1)/(k'*n)^2] * u'))
        % cvx_end
        %  
        %  
        %    Dr = subs(Dr,q, thr_d * rt_d)
    end
end